! This is a 1D advection example using square initial condition and periodic
! boundary condition for leap-frog finite difference scheme.
!
! Li Dong <dongli@lasg.iap.ac.cn>
!
! - 2018-03-20: Initial creation.
! - 2023-03-20: Add modified Robert-Asselin time filter (alpha=1 is original form, alpha=0 is no filter).
!               See: Williams, Paul D., 2009: A Proposed Modification to the Robert–Asselin Time Filter. Mon. Wea. Rev.
! - 2024-03-23: Use adv_1d_test_case_mod.

program leap_frog_adv_1d_case

  use adv_1d_test_case_mod

  implicit none

  real, allocatable :: rho(:,:)     ! Tracer density being advected at cell centers
  real, allocatable :: flx(:)       ! Flux at cell interfaces
  real d                            ! Time filter displacement
  real, parameter :: alpha = 1      ! Modified time filter weight
  integer, parameter :: ns = 1      ! Stencil width
  integer :: middle = 3, i
  character(30), parameter :: scheme = 'leap_frog'

  namelist /params/ nx, nt, dt, u

  call get_command_argument(1, namelist_path)
  inquire(file=namelist_path, exist=is_exist)
  if (is_exist) then
    open(10, file=namelist_path)
    read(10, nml=params)
    close(10)
  end if

  allocate(rho(1-ns:nx+ns,3))
  allocate(flx(1-ns:nx+ns  ))

  call adv_1d_test_case_init('square', ns, rho(:,old))
  rho(:,middle) = rho(:,old)
  call output(scheme, 0, ns, nx, x, rho(:,old))

  ! Run integration.
  print *, time_step, sum(rho(1:nx,old))
  do while (time_step < nt)
    call leap_frog(rho(:,middle), flx)
    do i = 1, nx
      rho(i,new) = rho(i,old) - dt / dx * (flx(i) - flx(i-1))
      ! Do time filtering.
      d = 0.1 * (rho(i,old) - 2 * rho(i,middle) + rho(i,new))
      rho(i,middle) = rho(i,middle) + alpha * d
      rho(i,new) = rho(i,new) + (alpha - 1) * d
    end do
    call apply_bc(ns, nx, rho(:,new))
    ! Change time indices.
    i = old; old = middle; middle = new; new = i
    time_step = time_step + 1
    call output(scheme, time_step, ns, nx, x, rho(:,old))
    print *, time_step, sum(rho(1:nx,old))
  end do

  call adv_1d_test_case_final()

  deallocate(rho)
  deallocate(flx)

contains

  subroutine leap_frog(q, f)

    real, intent(in ) :: q(1-ns:nx+ns)
    real, intent(out) :: f(1-ns:nx+ns)

    integer i

    do i = 1, nx
      f(i) = u * (q(i+1) + q(i))
    end do
    call apply_bc(ns, nx, f)

  end subroutine leap_frog

end program leap_frog_adv_1d_case
